% www.sherrytowers.com/mle_introduction.pdf
% f is the likelihood function

function [H_small, H, H_big, G] = own_hessian_v2(f, x0)

% Basic idea:
% - Log likelihood
% Partial derivatives
% then look at second derivatives of log likelihood
% Fisher information is negative of the hessian
% Variance covariance matrix is the inverse of the Fisher information

    n = length(x0);
    H = zeros(n,n);
    H_small = zeros(n,n);
    H_big = zeros(n,n);
    
    f0 = f(x0);
    G = own_gradient(f, x0, f0);
    G

    
    for i=1:n
        [i n]
        for j=1:i
            
            %[i j]
            
            % set the offsets for each of the three computations
            offset = abs(0.01 * x0(j)); % 1% of the magnitude of the estimate
            i_offset = abs(0.01 * x0(i)); % for the inner part of cross partial calculations
            
            offset_small = abs(0.005 * x0(j)); % 1% of the magnitude of the estimate
            i_offset_small = abs(0.005 * x0(i)); % for the inner part of cross partial calculations
            
            offset_big = abs(0.2 * x0(j)); % 1% of the magnitude of the estimate
            i_offset_big = abs(0.2 * x0(i)); % for the inner part of cross partial calculations
            
            if i==j  
                first_offset = zeros(size(x0));
                first_offset(i) = offset;
                H(i,i) = (f(x0+first_offset)-2*f0+f(x0-first_offset))/offset^2;                
                
                first_offset = zeros(size(x0));
                first_offset(i) = offset_small;
                H_small(i,i) = (f(x0+first_offset)-2*f0+f(x0-first_offset))/offset_small^2;                
                
                first_offset = zeros(size(x0));
                first_offset(i) = offset_big;
                H_big(i,i) = (f(x0+first_offset)-2*f0+f(x0-first_offset))/offset_big^2;  
                
            else % cross partial
                
                second_offset = zeros(size(x0));
                second_offset(j) = offset;
                H(i,j) = (own_deriv(f, x0 + second_offset, i, i_offset) -...
                             G(i))/ offset;

                second_offset = zeros(size(x0));
                second_offset(j) = offset_small;
                H_small(i,j) = (own_deriv(f, x0 + second_offset, i, i_offset_small) -...
                             G(i))/ offset_small;
                         
                second_offset = zeros(size(x0));
                second_offset(j) = offset_big;
                H_big(i,j) = (own_deriv(f, x0 + second_offset, i, i_offset_big) -...
                             G(i))/ offset_big;
          
            end       
        end
    end 
       
    
end